#ifndef MESH_H
#define MESH_H
#include <map>
#if defined(_MSC_VER)
	#include <hash_map>
#elif defined(__GNUC__)
	#include <ext/hash_map>
#endif

#include <fstream>
#include <iostream>
#include <sstream>
#include <string>

class HalfEdge ;
class Edge;
class Face;
class Vertex; 

class HalfEdge 
{
public:
	//HalfEdge (): next(NULL), prev(NULL), flip(NULL), face(NULL), from_v(NULL), edge(NULL){}
	HalfEdge * next, * prev, * flip;
	Face * face;
	Vertex * from_v ;
	Edge * edge; 
};

/*
class SamplePoint
{
public:
	int_RGB color ;
};
*/
class Vertex
{
public:
	Vertex ():locked(false), idx(-1) {}
	Vertex(
        const Point3 & p, const Point2 & t
    ) : point(p), tex_coord(t), locked(false), idx(-1) { }
public:
	class VFIter 
	{
	public:
		VFIter (Vertex * v)
		{
			he = v->he; 
		};
		VFIter (HalfEdge * edge)
		{
			he = edge; 
			begin = edge; 
		}
		VFIter & operator ++()
		{
			if (this->he->flip == NULL)	
				this->he = begin ;
			else
				this->he = he->flip->next ;
			return (*this) ;
		}
		Face * operator * ()
		{
			return (this->he->face);
		}
		bool operator ==(const VFIter & other)
		{
			return(he == other.he);
		}

		bool operator !=(const VFIter & other)
		{
			return(he != other.he);
		}

	private :
		HalfEdge * he; 
		HalfEdge * begin ;
	};
	//This two iterator will deal with the mesh has boundary
	VFIter vfbegin_safe()
	{
		//Check wheter the vertex is on boundary or not
		HalfEdge * iter  = this->he;
		HalfEdge * iterEnd  = this->he;
		do
		{
			if (iter->prev->flip == NULL)
			{
				bound = iter ;
				return VFIter (bound) ;
			}
			iter = iter->prev->flip ;
		}while (iter != iterEnd) ;
		bound = iterEnd ;
		return VFIter (bound) ;
	}
	VFIter vfend_safe()
	{
		return VFIter (bound);
	}
	VFIter vfbegin ()
	{
		return VFIter (this) ;
	}
	VFIter vfend ()
	{
		return VFIter (this->he->prev->flip);
	}
public:
	class VVIter 
	{
	public:
		VVIter (Vertex * v)
		{
			he = v->he; 
		};
		VVIter (HalfEdge * edge)
		{
			he = edge; 
		}
		VVIter & operator ++()
		{
			this->he = he->flip->next ;
			return (*this) ;
		}
		Vertex * operator*()
		{
			return (this->he->flip->from_v);
		}
		bool operator ==(const VVIter & other)
		{
			return(he == other.he);
		}

		bool operator!=(const VVIter & other)
		{
			return(he != other.he);
		}
		const  HalfEdge * get_he () const {return this->he ;}
	private :
		HalfEdge * he;  
	};
public:
	VVIter vvbegin ()
	{
		return VVIter (this) ;
	}
	VVIter vvend ()
	{
		return VVIter (this);
	}
public:
	const Point3 & normal () const 
	{
		return this->normal_vector ;
	}

public: 
	Point3 point; 
	Point2 tex_coord; 
	Point3 normal_vector ;
	HalfEdge * he; 
	unsigned idx ;

	unsigned num_neigh ;
	Point sum_neigh ;
	float qua_sum_neigh ; 

	//Used for computing conformal mapping
	bool locked ;

	//Used for mesh with boundary
	HalfEdge * bound ;
};

class Layer
{
public:
	Layer () : resolution (64) 
	{
		pixel.resize ( (resolution +1)  *(resolution+2)/2);
	}
	Layer (unsigned short r): resolution (r) 
	{
		pixel.resize ( (resolution +1)  *(resolution+2)/2);
	}
	Layer (const int_RGB & default_color) : resolution (16)
	{
		pixel.resize ( (resolution +1)  *(resolution+2)/2);
		for (int i = 0; i < (resolution +1)  *(resolution+2)/2; ++i)
		{
			pixel[i]= default_color ;
		}
	}
	void upSampling () 
	{
		std::vector<int_RGB> old_pixel (pixel);
		unsigned short res = this->resolution ;
		res *= 2 ;
		this->pixel.resize( (res+1) * (res+2) / 2);
		this->resolution = res; 

		for (unsigned i = 0; i <=res; ++i)
		{
			for (unsigned j = 0; j <= res - i ; ++j)
			{
				int k = res - i - j;
				if (i % 2 == 0 && j % 2 == 0 && k % 2 == 0)
				{
					pixel[barycentric2index(res, i,j)] =
						old_pixel[barycentric2index(res/2, i/2,j/2)] ;
					continue;
				}
				//Fix i, two neighbors are (i,j-1,k+1) and (i, j+1,k-1)
				if (i%2 == 0)
				{
					//pixel[barycentric2index(res, i, j)] =
					//	old_pixel[barycentric2index(res/2, i/2,(j-1)/2)] * 0.5f 
					//	+ old_pixel[barycentric2index(res/2, i/2,(j+1)/2)] * 0.5f;
					pixel[barycentric2index(res, i, j)] =
					linear_blend(old_pixel[barycentric2index(res/2, i/2,(j-1)/2)] ,
					old_pixel[barycentric2index(res/2, i/2,(j+1)/2)]);
					continue;
				}
				//Fix j, two neighbors are (i-1,j,k+1) and (i+1, j,k-1)
				if (j%2 == 0)
				{
					//pixel[barycentric2index(res, i, j)] =
					//	old_pixel[barycentric2index(res/2, (i-1)/2,(j)/2)] * 0.5f 
					//	+ old_pixel[barycentric2index(res/2, (i+1)/2,(j)/2)] * 0.5f;

					pixel[barycentric2index(res, i, j)] =
					linear_blend(	old_pixel[barycentric2index(res/2, (i-1)/2,(j)/2)]
						, old_pixel[barycentric2index(res/2, (i+1)/2,(j)/2)] );
					continue;
				}
				//Fix k, two neighbors are (i-1,j+1,k) and (i+1, j-1,k)
				if (k%2 == 0)
				{
					//pixel[barycentric2index(res, i, j)] =
					//	old_pixel[barycentric2index(res/2, (i+1)/2,(j-1)/2)] * 0.5f 
					//	+ old_pixel[barycentric2index(res/2, (i-1)/2,(j+1)/2)] * 0.5f;
					pixel[barycentric2index(res, i, j)] =
						linear_blend( old_pixel[barycentric2index(res/2, (i+1)/2,(j-1)/2)] , 
						old_pixel[barycentric2index(res/2, (i-1)/2,(j+1)/2)] );
					
					continue;
				}
			}
		}
	}
	void downSampling () 
	{
		unsigned short res = this->resolution ;
		res /= 2;
		for (unsigned i = 0 ; i <= res; ++i)
		{
			for (unsigned j = 0; j <= res - i; ++j)
			{
				this->pixel[barycentric2index(res, i,j)] = 
					this->pixel [barycentric2index (res * 2, i *2, j * 2)] ;
			}
		}
		this->resolution = res; 
		this->pixel.resize((res+1) * (res+2) / 2);
	}
	inline int_RGB & operator [](unsigned i)  
	{
		return pixel [i] ;
	}

	inline const unsigned short getResolution() const 
	{
		return this->resolution ;
	}
protected:
	unsigned short resolution; 
	std::vector<int_RGB> pixel ;
};

class Face 
{
public:
	typedef Texture::frgba_type frgba_type ;
	typedef math::vector<float, 3> point_type ;
	Face ():  selected (false), covered(false), reached (false)
	{
		Layer layer;
		// We create three layers. Here. 
		// The first layer is a white background.
		// The second and third layer is transparent. 
		Layer white_layer(int_RGB (255,255,255,255));

		//this->m_layers.push_back(white_layer);
		this->m_layers.push_back(layer);
		this->m_layers.push_back(layer);
		this->m_layers.push_back(layer);

		this->m_activeLayer = 1; 
		this->blenderFlag = false;
		this->blendLayers(); 

		strongEdge[0] = strongEdge[1] = strongEdge[2] = false;
	}

	inline void  Face2VertexIdx ( unsigned idx[3])
	{
		idx[0] = this->he->from_v->idx ;
		idx[1] = this->he->next->from_v->idx ;
		idx[2] = this->he->next->next->from_v->idx ;
	}
	inline void  Face2TexCoord (Point2 v[3])
	{
		v[0] = this->he->from_v->tex_coord ;
		v[1] = this->he->next->from_v->tex_coord ;
		v[2] = this->he->next->next->from_v->tex_coord ;
	}
	inline void  Face2VertexCoord (Point3 v[3])
	{
		v[0] = this->he->from_v->point ;
		v[1] = this->he->next->from_v->point ;
		v[2] = this->he->next->next->from_v->point ;
	}
	inline void  Face2Vertex (Vertex * v[3])
	{
		v[0] = this->he->from_v ;
		v[1] = this->he->next->from_v;
		v[2] = this->he->next->next->from_v;
	}
	inline const Vector3 & normal () const  
	{
		return this->normalVector ; 
	}
	void compute_param()
	{
		Point p[3];
		Face2VertexCoord(p);
		float nk, nh, nc;
		compute_tex_param(p, nk, nh, nc);
		mk = nk ;
		mh = nh ;
		mc = nc ;
	}
	inline void compute_tex_param (Point p[3], float & k, float & h, float &c)
	{
		//compute c & k ;
		float a = (p[1]-p[2]).length(); 
		float b = (p[2]-p[0]).length();
		c = (p[0]-p[1]).length();
		k = (b*b + c * c - a * a) / (2.0f * c) ; 

		//compute h ;
		float cos_a = (b*b + c* c - a* a)/ (2.0f * b * c) ;
		float sin_a = sqrt (1.0f - cos_a * cos_a) ;
		h = b * sin_a ;

		//computer c ;

		if (h >= c) {k /= h ; c /= h; h = 1.0f;}
		else {k /=c; h/= c; c = 1.0f;}
	}
	~Face() {}


	unsigned halfEdgeCount (HalfEdge * e)
	{
		HalfEdge * result = this->he ;
		for (unsigned i = 0; i < 3; ++i)
		{
			if (e == result) return i ;
			result = result->next ;
		}
		return 3; 
	}
	HalfEdge * halfEdge (unsigned count)
	{
		HalfEdge * result = this->he; 
		for (unsigned i = 0; i < count ; ++i) 
			result = result->next ;

		return result ;
	}
public: 
	//This group of function are used for painting .

	void setColor (const int_RGB & color, unsigned sample_idx)
	{
		//this->samples [sample_idx] = color ;
		this->m_layers[m_activeLayer][sample_idx] = color ;
	}
	const int_RGB & getColor (unsigned sample_idx)  
	{ return this->m_layers[m_activeLayer][sample_idx];}
	void upSampling () 
	{
		if (this->resolution  >= 64) return ; 

		for (unsigned i = 0 ;i < m_layers.size (); ++i)
		{
			m_layers[i].upSampling(); 
		}
		blendLayers();
	}
	void downSampling() 
	{
		if (this->resolution <= 1) return ;
		for (unsigned i = 0 ;i < m_layers.size (); ++i)
		{
			m_layers[i].downSampling();
		}
		blendLayers();
	}
	void setBlenderFlag (bool f = true)
	{
		this->blenderFlag = f; 
	}
	void setActiveLayer (unsigned short idx_layer) {this->m_activeLayer = idx_layer ;}
	const short getActiveLayer () const { return this->m_activeLayer;}

	void blendLayers ()
	{

		this->resolution = m_layers[0].getResolution();
		this->samples.resize ( (resolution + 1) * (resolution + 2) / 2);
		if (this->m_layers.size())
		{
			if (this->blenderFlag == false)
			{
				for (unsigned i = 0; i < this->samples.size(); ++i)
				{
					int_RGB color (255,255,255,255);
					if ( this->m_layers[m_activeLayer][i] == int_RGB(-1,-1,-1, 255))
					{
						this->samples[i] = int_RGB (255,255,255);
						continue ;
					}
					
					//Assume that all the layer has white background 
					this->samples[i] = solid_color (alpha_blend(m_layers[m_activeLayer][i],color));
						//solid_color (m_layers[m_activeLayer][i]) +
						//color * (1- m_layers[m_activeLayer][i].a /255.0f);
				}
			}
			else
			{
				for (unsigned i = 0; i < this->samples.size(); ++i)
				{
					bool flag = false;  
					int_RGB color (255,255,255,255);
					//for (int j = m_layers.size() - 1; j >= 0; --j)
					for (unsigned j = 0; j < m_layers.size(); ++j)
					{
						if (m_layers[j][i] == int_RGB (-1,-1,-1)) continue; 
						flag = true; 
						color = solid_color (m_layers[j][i] ) +  color * (1 - m_layers[j][i].a /255.0f );

					}
					if (flag == false) color = int_RGB (255,255,255,255);
					this->samples[i] = color ;
				}
			}
		}
		return ;
	}

	// New CODE
	private :
		float k , h  , c;
		GLfloat compute_k(point_type p[3])
		{
			float a = distance(p[1], p[2]) ;
			float b = distance(p[2], p[0]) ;
			float c = distance(p[0], p[1]) ;
			float k = (b*b + c * c - a * a) / (2.0 * c) ; 
			return GLfloat (k); 
		}
		GLfloat compute_h(point_type p[3])
		{
			float a = distance(p[1], p[2]) ;
			float b = distance(p[2], p[0]) ;
			float c = distance(p[0], p[1]) ;
			float cos_a = (b*b + c* c - a* a)/ (2.0 * b * c) ;
			float sin_a = sqrt (1.0 - cos_a * cos_a) ;
			return (b * sin_a) ;
		}
		GLfloat compute_c(point_type p[3])
		{
			return distance(p[0] , p[1]) ;
		}
		void compute_normalized_tex_param (point_type point[3])
		{
			k = compute_k(point);
			h = compute_h(point);
			c = compute_c(point);
			if (h >=c) { k /= h ; c/=h ; h = 1.0f ; }
			else{ k /= c; h /= c; c = 1.0; }
			if (k < 0)
			{
				float ratio = c - k ;
				k /= ratio ; c /= ratio; h /= ratio ; 
			}
			if (k > 1.0)
			{
				float ratio = c + k ;
				k /= ratio ; c /=ratio ; h /=ratio ;
			}
			k *= 0.8; c *=0.8; h *=0.8; 
		}
	public:
		void write_tex_param ()
		{
			unsigned row = m_idx ;
			unsigned col = 0; 
			TexCoord texcoord ; 
			texcoord.x = col + 0.1f;  // Padding
			texcoord.y = row + 0.1f;  // Padding
			ptm->update (row, 0, frgba_type (this->resolution, 0.0f, 0.0f, 0.0f));

			if (k>=0)
			{
				ptm->update(row, 1, frgba_type (texcoord.x, texcoord.y, 0.0f, 0.0f));
				ptm->update (row, 2, frgba_type (texcoord.x + c, texcoord.y , 0.0f, 0.0f));
				ptm->update (row, 3, frgba_type (texcoord.x + k, texcoord.y + h, 0.0f, 0.0f));
			}
			else if (k < 0)
			{
				ptm->update(row, 1, frgba_type (texcoord.x - k, texcoord.y, 0.0f, 0.0f)); 
				ptm->update(row, 2, frgba_type (texcoord.x - k + c, texcoord.y, 0.0f, 0.0f )); 
				ptm->update(row, 3, frgba_type (texcoord.x, texcoord.y + h, 0.0f, 0.0f)); 

			}

			for (unsigned i = 4; i < (this->resolution + 1 ) * (this->resolution + 2) / 2 + 4; ++i ) 
			{
				//                              ptm->update (row, i, frgba_type (1.0f, 1.0f, 1.0f, 1.0f));

				if (i % 3 == 1)
					ptm->update (row, i, frgba_type (1.0f, 0.0f, 0.0f, 1.0f));
				else if (i %3 == 2)
					ptm->update (row, i, frgba_type (0.0f, 1.0f, 0.0f, 1.0f));
				else 
					ptm->update (row, i, frgba_type (0.0f, 0.0f, 1.0f, 1.0f));
			}
		}
public:
	void draw(point_type point[3])
	{
		unsigned row = m_idx ;
		unsigned col = 0; 
		TexCoord texcoord ; 
		texcoord.x = col + 0.1f;  // Padding
		texcoord.y = row + 0.1f;  // Padding 

		glBegin(GL_TRIANGLES);
		point_type normal (normalVector.x, normalVector.y, normalVector.z) ;
		glNormal3f (normal.x, normal.y , normal.z); 
		if (k >=0 )
		{
			glTexCoord2f(texcoord.x,                texcoord.y); glVertex3f(point[0].x, point[0].y, point[0].z);
			glTexCoord2f(texcoord.x + c ,   texcoord.y ); glVertex3f(point[1].x, point[1].y, point[1].z);
			glTexCoord2f(texcoord.x + k ,   texcoord.y + h); glVertex3f(point[2].x, point[2].y, point[2].z);
		}
		else
		{
			glTexCoord2f(texcoord.x - k,            texcoord.y); glVertex3f(point[0].x, point[0].y, point[0].z);
			glTexCoord2f(texcoord.x - k + c ,       texcoord.y ); glVertex3f(point[1].x, point[1].y, point[1].z);
			glTexCoord2f(texcoord.x ,       texcoord.y + h); glVertex3f(point[2].x, point[2].y, point[2].z);
		}
		glEnd();        
	}
	void set_texture (Texture * tm) {this->ptm = tm ;}
	// new code end
public: 
	unsigned m_idx ;
	Texture * ptm ;
	HalfEdge * he;
	unsigned resolution; 
	int_RGB color; 
	bool selected;
	bool covered ; // For brush covered on the mesh 
	bool reached ; // For BFS on the mesh 
	

	//The data samples are merged from different layers
	std::vector<int_RGB> samples ; 

	//Layers supported
	std::vector<Layer> m_layers ;
	short m_activeLayer; 
	bool blenderFlag ;

	//Select edge support; 
	bool strongEdge[3];

	Vector3 normalVector ;


	float mk , tk;
	float mh, th; 
	float mc, tc;
} ;



class Edge 
{
public :
	HalfEdge * he ; 
	friend std::ostream & operator << (std::ostream & ostr, Edge * e)
	{
		ostr <<'('<<e->he->from_v->idx<<','; 
		ostr <<e->he->flip->from_v->idx<<')';
		return ostr;
	}
	float cost; 
};

typedef std::pair<int, int> HashKey;
typedef HashKey HalfEdgeHashKey ;
typedef HashKey EdgeHashKey; 

class FaceHashKey
{
public:
	FaceHashKey(unsigned a, unsigned b, unsigned c)
	{
		idx1 = a; idx2 = b; idx3 = c;
	}
public:
	unsigned idx1, idx2, idx3 ;
};


#if defined(_MSC_VER)


template<>
struct stdext::hash_compare<HashKey > {
    // the mean number of elements per bucket
    static const size_t bucket_size = 4;
    // the minimum number of buckets
    static const size_t min_buckets = 65535*16;
    // hash function
    size_t operator()(const HashKey & key) const {
        return (key.first * key.second);
    }
    // compare function
    bool operator()(const HashKey & lhs, const HashKey & rhs) const {
        return lhs < rhs;
    }
};

template<>
struct stdext::hash_compare<FaceHashKey > {
    // the mean number of elements per bucket
    static const size_t bucket_size = 4;
    // the minimum number of buckets
    static const size_t min_buckets = 65535*2;
    // hash function
    size_t operator()(const FaceHashKey & key) const {
		return key.idx1 * key.idx2 + key.idx3 ;
    }
    // compare function
    bool operator()(const FaceHashKey & lhs, const FaceHashKey & rhs) const {
		return (lhs.idx1 < rhs.idx1 || 
			(lhs.idx1 == rhs.idx1 && lhs.idx2 < rhs.idx2) ||
			(lhs.idx1 == rhs.idx1 && lhs.idx2 == rhs.idx2 && lhs.idx3 < rhs.idx3)) ;
    }
};


#elif defined(__GNUC__)
	struct EdgeHashOperator {
		size_t operator()(const EdgeHashKey & key) const {
        return (key.first * key.second);
    }
	};
	struct EdgeHashEqual {
		bool operator()(const EdgeHashKey & lhs, const EdgeHashKey & rhs) const {
			return lhs == rhs;
    }
	};

	struct FaceHashOperator {
		size_t operator()(const FaceHashKey & key) const {
        return (key.idx1 * key.idx2 + key.idx3);
    }
	};
	struct FaceHashEqual {
		bool operator()(const FaceHashKey & lhs, const FaceHashKey & rhs) const {
			return ( lhs.idx1 == lhs.idx1 
				  && lhs.idx2 == lhs.idx2 
				  && lhs.idx3 == lhs.idx3 );
    }
};

#endif
inline HalfEdgeHashKey halfEdgeHashKey (HalfEdge * he)
{
	unsigned id1 = he->from_v->idx ;
	unsigned id2 = he->next->from_v->idx; 
	return std::make_pair (id1, id2);
}
inline EdgeHashKey edgeHashKey (unsigned id1, unsigned id2)
{
	return (id1 < id2 ? std::make_pair(id1, id2) : std::make_pair(id2, id1)) ;
}
inline EdgeHashKey edgeHashKey(Edge * edge)
{
	if (edge == NULL) 
	{
		std::cout<<"INVALID EDGE"<<std::endl;
	}
	unsigned id1 = edge->he->from_v->idx;
	unsigned id2 = edge->he->next->from_v->idx ;
	return edgeHashKey(id1, id2);
}
inline FaceHashKey faceHashKey (unsigned id1, unsigned id2, unsigned id3)
{
	if (id1 < id2  && id1 < id3) return FaceHashKey(id1, id2, id3 );
	if (id2 < id1  && id2 < id3) return FaceHashKey(id2, id3, id1 );
	//if (id3 < id1  && id3 < id2) 
	return FaceHashKey(id3, id1, id2 );
}

inline FaceHashKey faceHashKey (Face * face)
{
	unsigned id1 = face->he->from_v->idx; 
	unsigned id2 = face->he->next->from_v->idx ;
	unsigned id3 = face->he->next->next->from_v->idx ;

	return faceHashKey (id1, id2, id3 );
}
inline void adjust_vertex (Vertex * v1, Vertex * v2)
{
	//std::vector <Vertex * > vertices; 
	v1->point = Point(v1->sum_neigh + v2->sum_neigh) / float(v1->num_neigh + v2->num_neigh); 
	v1->num_neigh +=  v2->num_neigh;
	v1->qua_sum_neigh += v2->qua_sum_neigh;
	v1->sum_neigh = v1->sum_neigh + v2->sum_neigh ;
}

inline bool is_edge_safe (Edge * edge)
{
	Vertex * v1 = edge->he->from_v ;
	Vertex * v2 = edge->he->next->from_v; 
	unsigned count = 0;
	Vertex::VVIter iter1 = v1->vvbegin ();
	Vertex::VVIter end1 = v1->vvend();
	do 
	{
		Vertex::VVIter iter2 = v2->vvbegin ();
		Vertex::VVIter end2 = v2->vvend();
		do
		{
			//std::cout<<"TEST"<<std::endl;
			if ( (*iter1)->idx == (*iter2)->idx ) count ++; 
			//std::cout<<(*iter2)->idx<<' ';
			++iter2;
		}while (iter2 != end2);
		++iter1;
	}while (iter1 != end1);
	return (2 == count);
}



#if defined(_MSC_VER)
	typedef stdext::hash_map <unsigned, Vertex * > VERTEX_CONTAINER ;
	typedef stdext::hash_map <EdgeHashKey, Edge * > EDGE_CONTAINER ;
	typedef stdext::hash_map <FaceHashKey, Face * > FACE_CONTAINER ;
	typedef stdext::hash_map <HashKey, HalfEdge * > HALFEDGE_CONTAINER ;
	
#elif defined(__GNUC__)
	typedef __gnu_cxx::hash_map <unsigned, Vertex * > VERTEX_CONTAINER ;
	typedef __gnu_cxx::hash_map <EdgeHashKey, Edge * , EdgeHashOperator, EdgeHashEqual> EDGE_CONTAINER ;
	typedef __gnu_cxx::hash_map <FaceHashKey, Face * , FaceHashOperator, FaceHashEqual> FACE_CONTAINER ;
	typedef __gnu_cxx::hash_map <HashKey, HalfEdge * , EdgeHashOperator, EdgeHashEqual> HALFEDGE_CONTAINER ;

#endif

	typedef VERTEX_CONTAINER::iterator VERTEX_ITERATOR ;
	typedef EDGE_CONTAINER::iterator EDGE_ITERATOR ;
	typedef FACE_CONTAINER::iterator FACE_ITERATOR ;
	typedef HALFEDGE_CONTAINER::iterator HALFEDGE_ITERATOR ;



class Mesh 
{
public: 
	friend std::istream & operator >> (std::istream & istr, Mesh & mesh)
	{
		std::string line;
		std::string tag; 
		std::string slash;
		Vertex * v; 
		unsigned count = 0; 
		unsigned tex_count = 0;

		while (std::getline(istr, line))
		{
			std::stringstream sstr(line.c_str());
			sstr >> tag;
			if (tag == "v")
			{	
				v = new Vertex;
				sstr >> v->point.x >> v->point.y >> v->point.z ;
				//v->point.z = - v->point.z;
				v->he = NULL; 
				v->idx = count ;
				mesh.vertices[count]=v;
				++ count ;
			}
			else if (tag == "vt")
			{
				Vector2 tex_coord;
				sstr>>tex_coord.x >> tex_coord.y ;
				mesh.vertices[tex_count]->tex_coord = tex_coord ;
				++ tex_count ;
			}
			else if (tag == "f")
			{
				std::vector <unsigned > v_id ;

                while(sstr) {
                    std::string s ;
                    sstr >> s ;
                    if(s.length() > 0) {
                        std::stringstream v_input(s.c_str()) ;
                        int index ;
                        v_input >> index ;
						v_id.push_back(index - 1) ;
                        char c ;
                        v_input >> c ;
                        if(c == '/') {
                            v_input >> index ;
                        }
                    }
                }
				mesh.construct(v_id[0], v_id[1], v_id[2]) ;
			}
		}
		for (VERTEX_CONTAINER::iterator iter = mesh.vertices.begin(); 
			iter!= mesh.vertices.end(); ++iter)
		{
			mesh.compute_normal((*iter).second); 
		}
		
		return istr;
	}
	friend std::ostream & operator << (std::ostream & ostr, Mesh & mesh)
	{
		for (unsigned i = 0 ;i < mesh.vertices.size(); ++i)
		{
			ostr <<"v " << mesh.vertices[i]->point<<std::endl;
		}
		for (unsigned i = 0 ;i < mesh.vertices.size(); ++i)
		{
			ostr <<"vt " << mesh.vertices[i]->tex_coord.x<<' '<< mesh.vertices[i]->tex_coord.y<<std::endl;
		}
		for (unsigned i = 0; i < mesh.v_faces.size(); ++i)
		{
			HalfEdge * he  = mesh.v_faces[i]->he; 
			ostr <<"f "<<he->from_v->idx + 1<<'/'<<he->from_v->idx + 1<<' '
				<<he->next->from_v->idx + 1<<'/'<<he->next->from_v->idx + 1<<' '
				<<he->next->next->from_v->idx + 1<<'/'<<he->next->next->from_v->idx + 1<<std::endl;
		}
		return ostr; 
	}
	void finalize ()
	{
		for (VERTEX_ITERATOR it = vertices.begin();
			it != vertices.end(); ++it)
		{
			if (it->second) delete (it->second);
			it->second = NULL; 			
		}
		std::cerr<<"Remove vertices"<<std::endl;

		for (FACE_ITERATOR it = faces.begin();
			it != faces.end(); ++it)
		{
			if (it->second) delete (it->second);
			it->second = NULL; 			
		}
		std::cerr<<"Remove faces"<<std::endl;
		
		for (EDGE_ITERATOR it = edges.begin(); 
			it !=edges.end(); ++it)
		{
			if (it->second) delete it->second; 
			it->second = NULL; 
		}
		std::cerr<<"Remove edges"<<std::endl;
		for (HALFEDGE_ITERATOR it = hash_he.begin(); 
			it != hash_he.end(); ++it)
		{
			if (it->second) delete it->second; 
			it->second = NULL; 
		}
		std::cerr<<"Remove Halfedges"<<std::endl;
	}
	~Mesh()
	{
		finalize();
	}


	void evaluate_edge_cost (Edge * edge)
	{
			Vertex * v1 = edge->he->from_v ;
			Vertex * v2 = edge->he->next->from_v ;

			Point p =  (v1->sum_neigh + v2->sum_neigh) / float(v1->num_neigh + v2->num_neigh) ;
			/*
			edge->cost = (v1->num_neigh + v2->num_neigh)*(p*p) - 
				p *(((v1->sum_neigh + v2->sum_neigh)) * 2.0) 
				+ v1->qua_sum_neigh + v2->qua_sum_neigh;
			*/
	}
	void evaluate_edge_around_vertex(Vertex * v)
	{
		HalfEdge * iter = v->he; 
		HalfEdge * end = iter;
		do
		{
			evaluate_edge_cost(iter->edge);
			iter = iter->flip->next ;
		} while(iter != end);
	}
	void prepare_simp_data()
	{
		for (VERTEX_ITERATOR v_iter = vertices.begin(); 
			v_iter != vertices.end(); ++v_iter)
		{
			unsigned count = 0;
			Point p(0.0, 0.0, 0.0);
			float sum = 0.0;
			Vertex::VVIter iter = (v_iter)->second->vvbegin();
			Vertex::VVIter end = (v_iter)->second->vvend();
			
			do 
			{
				++iter;
				++count;
				p = p + (*iter)->point; 
				sum += ((*iter)->point) * ((*iter)->point);
			} while (iter != end) ;
			v_iter->second->num_neigh = count;
			v_iter->second->sum_neigh = p; 
			v_iter->second->qua_sum_neigh = sum; 
		}
		for (EDGE_ITERATOR iter = edges.begin(); 
			iter != edges.end(); ++iter)
		{
			evaluate_edge_cost(iter->second);
		}
	}

private:
	void compute_normal (Vertex * v)
	{
		Point3 n (0.0f, 0.0f, 0.0f);
		Vertex::VFIter vfiter = v->vfbegin_safe();
		do 
		{
			n = n + (*vfiter)->normal();
			++ vfiter;
		} while (vfiter != v->vfend_safe());
		v->normal_vector = n.normalize(); 
	}

	void construct (const unsigned id1, const unsigned id2, const unsigned id3)
	{
		static unsigned idx_counter = 0; 
		
		Face * face = new Face; 
		HalfEdge * he12 = new HalfEdge ;
		HalfEdge * he23 = new HalfEdge ;
		HalfEdge * he31 = new HalfEdge ;
		he12->prev = he31 ; he12->next = he23; he12->face = face; he12->from_v = vertices[id1] ;
		he23->prev = he12 ; he23->next = he31; he23->face = face; he23->from_v = vertices[id2] ;
		he31->prev = he23 ; he31->next = he12; he31->face = face; he31->from_v = vertices[id3] ;

		hash_he[HashKey(id1 , id2)] = he12;
		hash_he[HashKey(id2 , id3)] = he23;
		hash_he[HashKey(id3 , id1)] = he31;
		
		vertices[id1]->he = he12; 
		vertices[id2]->he = he23; 
		vertices[id3]->he = he31; 		

		face->he =  he12 ;
		this->v_faces.push_back(face);

		faces[faceHashKey(id1, id2, id3)] = face;

		//Update the flips and Edges; 
		HALFEDGE_ITERATOR iter ;
		iter = hash_he.find(HashKey(id2, id1)) ;
		if (hash_he.end() == iter)
		{
			he12->flip = NULL;
			Edge * edge = new Edge;
			edge->he = he12 ;
			he12->edge = edge;
			if (id1 < id2) edges[HashKey(id1, id2)] = edge;
			else edges[HashKey(id2, id1)] = edge;
		}
		else
		{
			he12->flip = (*iter).second;
			(*iter).second->flip = he12 ;
			he12->edge = (*iter).second->edge ;
		}
		iter = hash_he.find(HashKey(id3, id2)) ;
		if (hash_he.end() == iter)
		{
			he23->flip = NULL;
			Edge * edge = new Edge;
			edge->he = he23 ;
			he23->edge = edge;
			if (id2 < id3) edges[HashKey(id2, id3)] = edge;
			else edges[HashKey(id3, id2)] = edge;
		}
		else
		{
			he23->flip = (*iter).second;
			(*iter).second->flip = he23 ;
			he23->edge = (*iter).second->edge ;
		}
		iter = hash_he.find(HashKey(id1, id3)) ;
		if (hash_he.end() == iter)
		{
			he31->flip = NULL;
			Edge * edge = new Edge;
			edge->he = he31 ;
			he31->edge = edge;
			if (id3 < id1) edges[HashKey(id3, id1)] = edge;
			else edges[HashKey(id1, id3)] = edge;
		}
		else
		{
			he31->flip = (*iter).second;
			(*iter).second->flip = he31 ;
			he31->edge = (*iter).second->edge ;
		}

		//update the normal
		Point3 p[3];
		face->Face2VertexCoord(p);
		Vector3 p1 = p[1] - p[0];
		Vector3 p2 = p[2] - p[0];
		face->normalVector = (p1^p2).normalize();
		face->compute_param();
		face->m_idx = idx_counter ;
		idx_counter++;
	}


public:
	void print_halfedges ()
	{
		std::cout<<"HALF EDGE SIZE :"<<edges.size()*2<<std::endl;
		for (EDGE_ITERATOR iter = edges.begin(); 
			iter != edges.end(); ++iter )
		{
			HalfEdge * he1 = iter->second->he ;
			HalfEdge * he2 = iter->second->he->flip; 
			if (he1->next->next->next != he1) std::cout<<"Invalid Halfedge"<<std::endl;
			if (he2->next->next->next != he2) std::cout<<"Invalid Halfedge"<<std::endl;
			std::cout<<'['<<he1->from_v->idx<<','<<he1->next->from_v->idx<<']'<<std::endl;
			std::cout<<'['<<he2->from_v->idx<<','<<he2->next->from_v->idx<<']'<<std::endl;
		}
	}
	void print_vertex_neighbor(Vertex * v)
	{
		std::cout<<std::endl<<"VERTEX NEIGHBOR : "<<v->idx<<std::endl;
		HalfEdge * iter = v->he; 
		HalfEdge * end = iter; 
		do
		{
			std::cout<<'['<<iter->from_v->idx <<' '<<iter->next->from_v->idx <<' '<<iter->next->next->from_v->idx <<']'<<std::endl;
			iter = iter->flip->next ;
		}while (iter != end);
		
	}
	void print_edges ()
	{
		std::cout<<"EDGE SIZE :"<<edges.size()<<std::endl;
		for (EDGE_ITERATOR iter = edges.begin(); 
			iter != edges.end(); ++iter )
		{
			
			std::cout<<'['<<iter->first.first<<','<<iter->first.second<<']'<<iter->second; 
			if (iter->first != edgeHashKey(iter->second->he->from_v->idx, iter->second->he->next->from_v->idx ))
			{
				std::cout<<" INVALID EDGE--";
				std::cout<<'['<<iter->first.first<<','<<iter->first.second<<']'<<
					iter->second->he->from_v->idx<<','<<iter->second->he->next->from_v->idx <<std::endl;
				exit(0);
			}
			std::cout<<std::endl;
		}
	}
	void print_vertex ()
	{
		for (VERTEX_ITERATOR iter = vertices.begin ();
			iter != vertices.end(); ++iter )
		{
			std::cout<<iter->second->idx << ": ("<< iter->second->num_neigh<<')';

			Vertex::VVIter iter_in = iter->second->vvbegin();
			Vertex::VVIter end = iter->second->vvend();
			do 
			{
				std::cout<<(*iter_in)->idx<<' ';
				++iter_in;
			} while (iter_in != end) ;
			std::cout<<std::endl;
			
		}
	}

	bool is_valid_mesh()
	{
		for (EDGE_ITERATOR iter = edges.begin(); 
			iter != edges.end(); ++iter )
		{
			HalfEdge * he1 = iter->second->he ;
			HalfEdge * he2 = iter->second->he->flip; 
			if (he1->next->next->next != he1)
			{
				std::cout<<"INVALID EDGE [CIRCLE1]:";
				return false;
			}
			if (he2->next->next->next != he2)
			{
				std::cout<<"INVALID EDGE [CIRCLE2]:";
				return false;
			}

			if (he1->next->from_v != he1->flip->from_v )
			{
				std::cout<<"\nINVALID EDGE [NEXT<->FLIP]1:"<<std::endl;
				print_vertex_neighbor(he1->from_v);
				print_vertex_neighbor(he1->next->from_v);
				print_vertex_neighbor(he1->flip->from_v);
				std::cout<<"NEXT "<<he1->next->from_v->idx<<' '<<he1->next->next->from_v->idx
					<<' '<<he1->next->next->next->from_v->idx<<std::endl;
				std::cout<<"FLIP "<<he1->flip->from_v->idx<<std::endl;
				return false;
			}
			if (he2->next->from_v != he2->flip->from_v )
			{
				std::cout<<"\nINVALID EDGE [NEXT<->FLIP]2:";
				print_vertex_neighbor(he2->from_v);
				print_vertex_neighbor(he2->next->from_v);
				print_vertex_neighbor(he2->flip->from_v);
				return false;
			}
			if (iter->first != edgeHashKey(iter->second->he->from_v->idx, iter->second->he->next->from_v->idx ))
			{
				std::cout<<std::endl;
				std::cout<<"INVALID EDGE [NEXT]:";
				std::cout<<'['<<iter->first.first<<','<<iter->first.second<<']'<<
					iter->second->he->from_v->idx<<','<<iter->second->he->next->from_v->idx <<std::endl;
				return false;
			}
			if (iter->first != edgeHashKey(iter->second->he->from_v->idx, iter->second->he->flip->from_v->idx ))
			{
				std::cout<<std::endl;
				std::cout<<" INVALID EDGE [FLIP]:";
				std::cout<<'['<<iter->first.first<<','<<iter->first.second<<']'<<
					iter->second->he->from_v->idx<<','<<iter->second->he->flip->from_v->idx <<std::endl;
				return false;
			}
		}
		return true;
	}
public:
	
	VERTEX_CONTAINER vertices;
	EDGE_CONTAINER edges; 
	FACE_CONTAINER  faces; 
	std::vector<Face * > v_faces;
	HALFEDGE_CONTAINER hash_he; 

};

class PseudoMesh 
{
public:
	PseudoMesh (Mesh * mesh): mesh (mesh)
	{

	}
public:
	inline unsigned pseudo2mesh_idx (unsigned pseudo_idx) const
	{
		return idx_pseudo2mesh[pseudo_idx];
	}
	inline unsigned mesh2pseudo_idx (unsigned mesh_idx) const
	{
		return idx_mesh2pseudo[mesh_idx];
	}

public: 
	inline Mesh * getMesh () {return this->mesh;}
	inline void  Face2MeshVertexIdx (Face * face, unsigned idx[3]) 
	{
		//mesh->Face2VertexIdx(face, idx );
		face->Face2VertexIdx(idx);
		idx[0] = pseudo2mesh_idx(idx[0]);
		idx[1] = pseudo2mesh_idx(idx[1]);
		idx[2] = pseudo2mesh_idx(idx[2]);
	}

	//The face pointers point to the original face
	std::vector<Face * > v_faces; 
	std::vector<unsigned> idx_mesh2pseudo ;
	std::vector<unsigned> idx_pseudo2mesh ;
private :
	Mesh * mesh; 



};
#endif
